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We present results from an extensive simulational study of 
the modulated liquid < > crystal transition in a 2-d charge- 
stabilized colloid subject to a 1-d laser field modulation com- 
mensurate with the crystalline phase. Contrary to some ear- 
lier simulational and experimental findings we do not find 
any reentrant liquid phase in our simulation. Furthermore 
the transition remains first-order (albeit weak, with a fairly 
large correlation length) even in the limit of infinite field, con- 
trary to mean-field predictions. In the modulated liquid phase, 
while the translational order decays exponentially, the bond 
orientational order is actually long ranged. 

PACS numbers : 64.70.Dv, 82.70.Dd, 05.70.Fh 



The freezing of charge-stabilized colloidal particles in 
the presence of a stationary modulation (laser) field has 
been of considerable research interest in recent years. Us- 
ing light scattering techniques Chowdhury, Ackerson and 
Clark |l[] showed that a charge stabilized colloidal liquid 
system, confined between two glass plates to form a sin- 
gle layer, when subjected to a commensurate stationary 
laser modulation (i.e. with the wavevector of the mod- 
ulating laser field tuned to be half the wavevector at 
which the direct correlation function S-'^\q) of the liq- 
uid develops its first peak in the absence of an exter- 
nal potential), freezes to form a triangular lattice with 
full two dimensional symmetry. They [Q also analyzed 
the phenomenon in terms of a simple Landau - Alexan- 
der - McTague (LAM) [p| theory and concluded that the 
transition from the modulated liquid to the crystal be- 
comes continuous (second-order) when the field strength 
is large. The continuous growth of intensities produced 
by all the density modes with increasing external field in- 
tensity was taken by Ackerson and Chowdhury |^ as an 
indication in favor of a second-order transition scenario. 
Later experimental studies involving direct microscopic 
observations Q and simulational studies using Monte- 
Carlo (MC) technique ||^ confirmed this phenomenon of 
laser induced freezing (LIFV 

Several authors Q, 0, [gl have studied this problem 
using density functional theory Q. The authors of ||] 
have shown from general symmetry grounds that for a 
suitably chosen modulation potential the free energy ex- 
pansion for the crystalline phase about the modulated 
liquid contains the relevant order parameters only in even 
powers. Hence there arises a possibility of change over 
of the first-order freezing transition at low external field 



strengths to a continuous transition for large enough ex- 
ternal field via a tricritical point. Their DFT results were 
in accordance with these symmetry arguments. 

In a later work Chakrabarti et al. |jlQ] studied a two 
dimensional colloidal system using MC simulation. The 
use of a standard procedure of looking at finite-size be- 
havior of the fourth cumulant of the energy |ll| seemed 
to confirm the existence of the tricritical point. Their 
study also found an intriguing reentrant modulated liq- 
uid phase - where increasing the external field strength 
actually melted the system instead of taking it towards 
the crystalline phase. Another unusual feature found in 
their simulation was that the order parameter jumped to 
zero at a lower value of the screening length compared to 
where the specific heat peak was seen. 

Very recently Wei et al. ||l2| reported an experimen- 
tal study of LIF and from the decay of pair correlation 
and the real space density plots concluded that their re- 
sults were in accordance with the results of reference [|0| , 
showing a reentrant liquid phase. 

In this letter we report the main results obtained from 
a vastly more extensive Monte-Carlo simulation than 
that in ref. |1^ of charge stabilized colloidal particles 
(with diameter 2R — 1.07 fim). They are assumed to 
be confined in a 2-d cell of size ^agL x a^L with pe- 
riodic boundary conditions and subjected to an exter- 
nal potential of the form U{r) = —VeCos{qox), with 
qq = 27r/(-^as), where Ug is the mean inter-particle sep- 
aration. The inter-particle interaction is modeled by the 
DLVO potential: 

(Ze)^ explnR) 2 exp{— Krij) 



U^,{r) 



1 + kR 



(1) 



Here Ze (Z — 7800) is the effective surface charge, e 
(=78) is the dielectric constant of the solvent and k is the 
inverse of the Debye screening length due to the counte- 
rions in the solvent. Details of the simulation will be 
published elsewhere ||l3|| . 

The phase diagram obtained by extrapolating the sim- 
ulation results to infinite system size as detailed below 
is shown in figure H^. It does not have the reentrance re- 
ported in pel ] . The transition is "strongly first-order" at 
zero field, signaled by a strong peak in the specific heat 
Cv [figure 0(a)] and sharp rise of the order parameters 
g , [figure |(c)], at K^a^ = 14.47. At larger field 
strengths, the peak in Cy becomes less prominent [figure 
0(d)] and the order parameter pd characterizing the mod- 
ulated liquid < — > crystal transition ||] , |l^ seems to go 
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FIG. 1. Phase diagram of laser induced freezing (extrapo- 
lated for infinite system size). 



to zero continuously [figure ||(f)], in apparent accordance 
with the mean field theory predictions 0] , [g[. Eow- 
ever, as we show below, in actuality the transition seems 
to remain first- order (albeit weak, with a fairly large cor- 
relation length) even in the limit of infinite field, contrary 
to the mean-field predictions. 

For a finite system, one expects that the transition 
will actually be rounded always. Hence one must look 
at the scaling of various quantities with large system size 
L to ascertain the order of the transition. In particu- 
lar the limiting value of the fourth cumulant of energy 
Vl = I — g'^ga'^a and its scaling behavior have been 
widely used to judge the order of the transition pi[ , 
ffo|. In case of a first-order transition, T4o = 2/3, the 



peak in Cy should diverge as L'^, while (k| 



3) and 



{Vl — Kx)) should scale as L~ . For small system sizes, 
there are corrections to this scaling behavior and we show 
below that at large field strengths these corrections be- 
come even more important. 

For a first-order transition scenario, assuming that the 
probability distribution for energy (per particle) is the 
sum of two Gaussian distributions centered around -E+ 
and _E_ corresponding to the disordered (modulated liq- 
uid) and the ordered (crystalline) phases, one can derive 
the following results, |11|, p^: 



CI = AL'^ + 

= K*+ ML- 



B + DL 



-d 



-2d 



C>(i-2d) 



(2) 
(3) 



where C£ refers to the height of the specific heat peak 
for a system of linear size L and A 



_ (B. 



E-) 



To 

leading order, the constants M, N and D are proportional 
to the inverse of ^§^, AF being the free energy difference 
between the liquid and the crystalline phases. 

Figures 0(a) and (c) show the scaling behavior of k^ 
(from the specific heat peak) for /3Ve = 0.0 and /3Ve = 2.0 
respectively. By fitting equation H (the solid line), we get 
K*^as = 14.34 ± 0.02(15.66 ± 0.03) for /3Ve = 0.0(2.0). 
C£ is plotted as function of L in figure ^(b) and (d). 
By fitting equation | we find l3{E+ - E-) ^ 0.071 ± 
0.0067(0.0044 ± 0.0031) for /3Ve = 0.0(2.0). 
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FIG. 2. Simulation results for 400 particles. The left hand 
panels (a-c) are for /3Ve = 0.0 and the right hand panels (d-f) 
are for (3Ve = 2.0. 



While the above analysis indicates that a first-order 
scenario is a likely possibility even at large /3Ve, the scal- 
ing analysis is done for results which span less than a 
decade in L (L = 8 to 30). Because the peak in Cy 
is rather broad and small in height, one needs to aver- 
age over several million MC steps to determine the peak 
position and peak height within acceptable errorbars. 
This limits the largest system we are able to simulate 
to L = 30. 
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FIG. 3. Finite size scaling behavior of the transition Kc(i) 
(a,b) and of the specific heat peak height (c,d) 

One can further argue that even though the transi- 
tion remains first order and does not show any reentrant 
melting till f3Ve — 2.0, there is always a possibility that 
the transition does become continuous at a larger field 



strength. To settle this we have carried out an even more 
extensive simulation at (3Ve = oo, so that the particles are 
confined to parallel lines. The resulting simplifications 
allow us to write a more efficient MC code, making sim- 
ulations upto L = 100 doable within our computing re- 
sources. The plot of K*^ as a function of L^ for /3Ve = oo, 
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FIG. 4. Finite size scaling for /3Vc = oo 

shown in Figure 0(a), has a maximum around L r^ SOa^. 
This non-monotonic behavior seem to suggest a large but 
finite length scale in the system. The existence of such a 
length scale is difficult to understand in the frame work 
of a second-order transition scenario, where, asymptoti- 
cally close to the phase transition, one should have just 
two length scales, the system size L (to which the "diverg- 
ing" correlation length ^ saturates) and the interaction 
scale (which is not much larger than as for the present 
case). But if the transition is weakly first-order with a 
large but finite f , one expects the finite-size scaling be- 
havior of the first-order transition discussed above to be 
valid only for L ^ S,- Figures Q(b) and ^(c) show that 
K^ and C2 for large L do follow a first-order like scal- 
ing [cf., eq.s (2) and (3)], with k*^ = 15.62 ± 0.05 and 
/3{E+ - £;_) ^ 0.0014 ± 0.0004; so does Vl which scales 
as L~^ as shown in figure 0(d). 

For /3Ve = oo we have computed the correlation 
functions for the coarse grained translational and bond- 
orientational order parameters, to be denoted by Griy) 
and Ge{y) respectively, for L — 70 and 100. In the mod- 
ulated liquid phase we obtain £_ by fitting Gt{v) to a 
decaying exponential. Figure p|(c) shows that S, actually 
becomes large even in the liquid phase as one goes to- 
wards K*^. But it saturates at around ISa^, indicating 
that the transition really is a weakly first-order transi- 
tion. 

Figures g(a) and (b) show G^iy) for the liquid and the 
solid phase near n*^. As can be seen from the figures, 
Ge{y) remains long ranged even in the liquid phase. A 
2-d liquid has a local hexagonal order. When one ap- 
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FIG. 5. correlation functions for /3Ve — oo 

plies a commensurate field, a unique direction (modulo 
60°) is picked up every where in the system. Thus the 
1-d modulation potential induces not only translational 
ordering, but also a long-range orientational order even 
in the modulated liquid phase. This is also understand- 
able from the general phcnomcnological mean-field free 
energy given in [n3, where a term linear in the bond- 
orientational order parameter arises as soon the transla- 
tional order parameter pf conjugate to Ve is turned on 
in the modulated liquid phase. Thus in the modulated 
liquid, Ge{y) goes to a constant value at large distances. 
What distinguishes the modulated liquid from the crys- 
talline phase is Griy), which decays exponentially in the 
modulated liquid phase but as a power law in the 2-d 
crystalline phase. The energy difference between the two 
phases becomes smaller as /3Ve is increased (though it 
never becomes strictly zero) and is responsible for the 
broad specific heat peak and the broad region over which 
the order parameter goes to zero. We believe that this 
probably led to the results reported in Q, where the 
true transition was missed, and a statistical fluctuation 
was mistakenly taken to be the Cy peak and hence the 
signature of a transition. Since different moments of en- 
ergy are not independent, a chance increase of Cy will 
be associated with a dip of Vl- Such artifacts arising 
from insufficient averaging are presumably absent in the 
present simulation due to the extensive averaging that 
has been carried out. This also explains why the earlier 
simulation ||lffl failed to see the order parameter fall to 
zero until much higher values of kGs ■ 

For moderate external fields, substantial diffusion pre- 
cedes the actual transition. While maintaining the trans- 
lational order over the full system, locally the system ac- 
commodates defects. In Figure g we show contour plots 
for the time averaged density as computed over different 
time spans from a Brownian dynamics simulation of 400 
particles for f]Ve = 2.0. The short-time averaged density 
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FIG. 6. Contour plots for the time averaged density from a 
Brownian dynamics simulation at /3Ve = 2.0 and kUs = 15.3. 
The 400 particle system was equilibrated in a run of 20 min- 
utes. Figures (a), (b) and (c) correspond to averaging over 
2 minute intervals starting from the th, 6 th and 16 th 
minute respectively after the system had equilibrated. Figure 
(d) corresponds to averaging over 20 minutes starting from 
equilibration. 



plots show near perfect crystalline order, with small re- 
gions where the particles have moved considerably from 
the lattice positions. If one averages for a long time (~ 20 
minutes of real time), the average density shows that the 
chance of finding a particle is nonzero everywhere along 
the lines. But we still find that the peaks of the average 
density maintain the crystalline order. Scattering exper- 
iments, which measure the structure at short time scales 
compared to the diffusion time, can differentiate this from 
the liquid phase via the existence of Bragg peaks. And in 
our simulation we can directly compute the order param- 
eter, which, for the above parameter values, is as large 
as 0.75! 

The claim in ||l^ of the observation of reentrant melt- 
ing into the modulated liquid phase at high (3Ve was 
based on the observed stronger decay of the envelope of 
g(r) than that theoretically predicted for the crystalline 
state |14|. The problem with this analysis is that the 
correct power-law decay in the envelope of g{r) sets in at 
much larger distances than for the order parameter corre- 
lation Gt (y) (which unfortunately is difficult to measure 
from experiments). In fig. H(d) we plot the envelope of 
g{y) and Gxiy) for f3Ve = oo, and kQs = 15.55 , i.e., in 
the crystalline phase but close to the transition . The en- 
velope of g{y) falls much more rapidly than Gt{jj) even 
upto y — 20 , and would yield a large exponent if fitted 
to a power law, in contrast to Gxir) which quickly at- 
tains a power law decay with exponent —0.3. We expect 
that for the same parameter values where |13| suggested a 
modulated liquid phase, light scattering experiments or 
order parameters calculated from Fourier transforming 



individual configurations prior to averaging will reveal 
the existence of a crystalline order. 

In conclusion we have presented new simulational evi- 
dence that the laser induced freezing in charge stabilized 
colloids remains first-order for an arbitrarily large laser 
field, in contrast to all theoretical predictions hitherto . 
The reentrant liquid phase suggested in jl^ and [Q is 
absent in the resulting phase diagram. Also we have 
shown that the modulated liquid phase has long range 
bond orientational order. The difference in the nature of 
transition between that predicted by the DFT and that 
found in our simulation is likely to be due to fluctuations, 
which are neglected in the DFT calculations in [6-8] ; the 
inclusion of such fluctuations is known to turn some tran- 
sitions flrst-order, eventhough they are predicted to be 
continuous by mean-field theory ||l^ ||l6[ . 
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